knitr::opts_knit$set(root.dir = '~/Documents/bgs_sims/')
knitr::opts_chunk$set(fig.width=16)

Libraries and functions

library(tidyverse)
library(ggridges)
library(viridis)
library(cowplot)
library(wesanderson)
orig_theme <- theme_set(theme_cowplot(font_size=20))

Data preparation functions

read_sim<-function(pidata,xidata,tajdata,roh=8.2e-10,neutral=FALSE,windows=20){
  #read data
  pi<-as.tibble(read.csv(pidata,header=T))
  print(nrow(pi))
  pi<-mutate(pi,stat=rep("pi",length(pi[,1])))
  xi<-as.tibble(read.csv(xidata,header=T)) 
  xi<-mutate(xi,stat=rep("xi",length(xi[,1])))
  tajD <- as.tibble(read.csv(tajdata,header=T))
  tajD<-mutate(tajD,stat=rep("tajD",length(tajD[,1])))
  #munge into single tibble of means
  data<-rbind(xi,pi,tajD) %>%
    `colnames<-`(c("gen", 1:windows, "sim", "stat")) %>%
    gather(window,value,-gen,-stat,-sim) %>% 
    group_by(gen,stat,window) %>%
    summarize(mean_val=mean(value,na.rm = T)) %>%
    mutate(centimorgen=(as.numeric(window)-1)*10000*roh*100 ,mean_val=mean_val)
  return(data)
}
#neutral=TRUE
#pidata <- "results/tennessen_neutral_pi.csv"
#xidata <-  "results/tennessen_neutral_xi.csv"
#tajdata <- "results/tennessen_neutral_tajD.csv"

Merge neutral and selected values

merge_sel_neutral <- function(sel_sim_data,neut_sim_data){
  merged <- merge(sel_sim_data,neut_sim_data,by=c('gen','stat','window','centimorgen'))
  merged <- mutate(merged,mean_val = if_else(stat=='tajD',mean_val.x-mean_val.y,mean_val.x/mean_val.y))
  names(merged) <- c('gen', 'stat', 'window', 'centimorgen','bgs','neutral','mean_val')
  merged
}

Plot functions

Plot selection landscape

simplot<-function(somedata,mystat,label){
  ggplot(filter(somedata,stat==mystat), aes(x=centimorgen, y=mean_val, group = gen,color=gen)) + 
  scale_color_viridis(name='Generation') +   
  background_grid(major = "xy", minor = "none") +
#  geom_smooth(se=F) + 
    geom_line() +
  ylab(label)  + xlab("Distance (cM)") #+
#  scale_x_continuous(breaks = seq(-1.5,1.5,0.5))
}

simplot_grid <- function(somedata,events,event_names,ratio=TRUE){
  if (ratio == TRUE){
    
    pi_lab <-  expression(bar(pi)/bar(pi)[neut])
    xi_lab <- expression(bar(xi)/bar(xi)[neut])
    d_lab <- expression(bar(D)-bar(D)[neut])
  }
  else{
    pi_lab <- expression(bar(pi))
    xi_lab <- expression(bar(xi))
    d_lab <- expression(bar(D))
  }
legend <- get_legend(simplot(somedata,"pi",expression(bar(pi)/bar(pi)[neut]))+
                       theme(legend.key.height= unit(2, "cm"),
                             legend.title = element_text(size=16),
                             legend.text = element_text(size=11)) +
                       scale_color_viridis(name='Generation',breaks = events, labels=event_names)
                     )
plot_grid(
plot_grid(
simplot(somedata,"pi",pi_lab)+theme(legend.position='none'),
simplot(somedata,"xi",xi_lab)+theme(legend.position="none"),
simplot(somedata,"tajD",d_lab)+theme(legend.position="none"),
ncol = 1,
align = 'hv'),
legend, rel_widths = c(85,15))
}

Function to plot ratios over time and by genetic distance

plot_ratios <- function(data,events,windows=20){
  window_list <- c(1,5,10,20)
comb <- data %>%
  filter(window%in%window_list,stat%in%c('pi','xi')) %>%
  mutate(centimorgen=abs(centimorgen)) %>%
  group_by(centimorgen,gen,stat) %>%
  summarise(mean_val=mean(mean_val))
start_val <- comb %>%
  filter(gen<0) %>%
  group_by(stat,centimorgen) %>%
  summarise(mean_val=mean(mean_val))

ggplot(comb,aes(gen,mean_val,color=stat)) +
  geom_smooth(method = 'loess',span=0.1)+
#  geom_line(size=1.5) +
  geom_hline(data=data.frame(yint=start_val$mean_val,centimorgen=start_val$centimorgen),aes(yintercept =yint))+
  labs(x=expression(paste('Generation (x ',N[anc],')')),y=expression(x/x[neut]),color='Statistics'
       #,caption ='Summary statistics at different distances (cM) from selected region.\nVertical lines mark demographic events. Horizontal lines denote value before demographic changes'
       ) + 
#  xlim(-0.1,1.25) +
  facet_grid(.~centimorgen) +
  geom_vline(xintercept = events, size=0.5, alpha=0.2) +
   scale_color_manual(values=c("#999999", "#E69F00"), 
                       breaks=c("pi", "xi"),
                       labels=c(expression(bar(pi)),expression(bar(xi)))) +
  theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
        strip.text.y = element_text(angle = 0,hjust = 0),
        plot.caption=element_text(size=12,hjust = 0)) +
#  scale_x_continuous(breaks=c(0,0.3,0.6,0.9,1.2)) +
NULL

}

Function to plot absolute values over time and by genetic distance

plot_over_time <- function(data,events,stats,windows=20){
  center <- 1
  window_list <- c(2,10,20)
  comb <- data %>%
  filter(window%in%window_list,stat==stats) %>%
  mutate(centimorgen=abs(centimorgen)) %>%
  group_by(centimorgen,gen,stat) %>%
  summarise(neutral=mean(neutral),bgs=mean(bgs)) %>%
  gather(selection,value,bgs,neutral) 

start_val <- comb %>%
  filter(gen<0) %>%
  group_by(stat,centimorgen,selection) %>%
  summarise(value=mean(value))
stats_label <- c(expression(bar(stats)))
ggplot(comb,aes(gen,value,color=selection)) +
  geom_line(size=1.5) +
  geom_hline(data=data.frame(yint=start_val$value,centimorgen=start_val$centimorgen,sel=start_val$selection),aes(yintercept =yint,color=sel))+
  labs(x=expression(paste('Generation (x ',N[anc],')')),y=stats,color='Selection')+
  facet_grid(.~centimorgen) +
  geom_vline(xintercept = events, size=0.5, alpha=0.2) +
   scale_color_manual(values=alpha(c("#ef8a62", "#67a9cf"),.9), 
                       breaks=c("neutral", "bgs"),
                       labels=c('Neutral','BGS')) +
  scale_x_continuous(breaks = seq(0,2,0.5))+
  theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
        strip.text.y = element_text(angle = 0,hjust = 0),
        plot.caption=element_text(size=12,hjust = 0))
}

Demographies real models

demog_torres <- read.csv('demographies/torres.csv')
demog_torres <- mutate(demog_torres,model='torres')
demog_ten <- read.csv('demographies/tennessen.csv')
demog_ten <- mutate(demog_ten,model='tennessen')
demog_maize <- read.csv('demographies/maize_018.csv')
demog_maize <- mutate(demog_maize,model='maize')
demog <- rbind(demog_torres,demog_ten,demog_maize)

# Torres events
anc_ext <- 0
OOA_torres <-   0.437017*2
euro_bneck_torres <- OOA_torres + 0.109481*2
torres_events <- c(anc_ext,OOA_torres,euro_bneck_torres)
torres_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck')
torres_event_sizes <- c(2.10709*18449,2.10709*18449,0.322295*18449)
torres_event_end<- c(1e5,1e5,2e3)

# Tennessen events
OOA_ten <- 0.264196*2
euro_bneck_ten <-OOA_ten+ 0.0763702*2
finel_ext <-euro_bneck_ten+ 0.0502841*2
tennessen_events <- c(anc_ext,OOA_ten,euro_bneck_ten,finel_ext)
tennessen_event_sizes <-  c(1.98*7310,1.98*7310,0.254609*7310,1.339863*7310)
tennessen_event_end<- c(1e5,1e5,1e5,1e3)
tennessen_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck','Final growth')

# maize events
maize_events <- c(0,0.12)
maize_event_end<- c(12278*0.04,12278*0.5)
maize_event_sizes <- c(12278,12278*2.98)
maize_event_names <- c('Domestication bottleneck','Maize present')

Extrapolate maize demog

x <- filter(demog_maize,generation>0)
fit <- lm(log(x$N)~(x$generation))
x$pred <- exp(fit$coefficients[1]+x$generation * fit$coefficients[2])
ggplot()+
  geom_line(data=x,aes(x=generation,y=pred),size=2) +
  geom_line(data=x,aes(x=generation,y=N),color='red')

exp(fit$coefficients[1]+0.18 * fit$coefficients[2])

Demography plot

ggplot(demog,aes(x=generation,y=N,color=model))+geom_line(size=2) + 
  scale_y_log10()+ 
  scale_color_discrete(breaks=c("tennessen", "torres","maize"),labels=c('Tennessen et al.','Torres et al.','Bessinger et al.'))+
  # annotate tennessen
  annotate("segment", x = tennessen_events, xend = tennessen_events, 
           y = tennessen_event_end, yend = tennessen_event_sizes, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = tennessen_events,  y = tennessen_event_end,
           label=tennessen_event_names,
           angle=c(90,90,90,0),hjust=0,vjust=c(0.5,0.5,0.5,1)) +
  # annotate torres
  annotate("segment", x = torres_events, xend = torres_events, 
           y = torres_event_end, yend = torres_event_sizes, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = torres_events,  y = torres_event_end,
           label=torres_event_names,
           angle=c(90,90,0),hjust=0,vjust=c(0.5,0.5,1)) +
  # annotate maize
   annotate("segment", x = maize_events, xend = maize_events, 
           y = maize_event_end, yend = maize_event_sizes, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = maize_events,  y = maize_event_end,
           label=maize_event_names,
           angle=0,hjust=0,vjust=1) +
  labs(x=expression(paste('Generation (x ',N[anc],')')),y= "Population size",color='Demography')

ggsave('figures/annotated_demographies.pdf',width = 16)
## Saving 16 x 5 in image

Europeans

Tennessen

Neutral data

tennessen_neutral <-read_sim("results/tennessen_neutral_pi.csv",
                             "results/tennessen_neutral_xi.csv",
                             "results/tennessen_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 610000
simplot_grid(tennessen_neutral,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

tennessen_bgs <-read_sim("results/tennessen_bgs_pi.csv",
                         "results/tennessen_bgs_xi.csv",
                         "results/tennessen_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 610000
simplot_grid(tennessen_bgs,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

tennessen <- merge_sel_neutral(tennessen_bgs,tennessen_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(tennessen,tennessen_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(tennessen,tennessen_events,'pi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'xi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/origVal_distance_summary_tennessen.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(tennessen,tennessen_events,tennessen_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/region_summary_tennessen.pdf',width = 12,heigh=9)

plot_ratios(tennessen,tennessen_events)

ggsave('figures/smooth_distance_summary_tennessen.pdf',width = 16,height = 7)

Torres

Neutral data

torres_neutral <-read_sim("results/torres_neutral_pi.csv",
                             "results/torres_neutral_xi.csv",
                             "results/torres_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2230000
simplot_grid(torres_neutral,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

torres_bgs <-read_sim("results/torres_bgs_pi.csv",
                         "results/torres_bgs_xi.csv",
                         "results/torres_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2230000
simplot_grid(torres_bgs,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

torres <- merge_sel_neutral(torres_bgs,torres_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(torres,torres_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(torres,torres_events,'pi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'xi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/origVal_distance_summary_torres.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(torres,torres_events,torres_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/region_summary_torres.pdf',width = 12,heigh=9)

plot_ratios(torres,torres_events)

ggsave('figures/smooth_distance_summary_torres.pdf',width = 16,height=7)

Maize

Neutral data

maize_neutral <-read_sim("results/maize_neutral_pi.csv",
                             "results/maize_neutral_xi.csv",
                             "results/maize_neutral_tajD.csv",roh = 8.2e-9,
                             neutral=TRUE)
## [1] 165000
simplot_grid(maize_neutral,maize_events,maize_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

maize_bgs <-read_sim("results/maize_bgs_pi.csv",
                         "results/maize_bgs_xi.csv",
                         "results/maize_bgs_tajD.csv",roh = 8.2e-9,
                         neutral=FALSE)
## [1] 165000
simplot_grid(maize_bgs,maize_events,maize_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

maize <- merge_sel_neutral(maize_bgs,maize_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(maize,maize_events,'pi',windows = 21))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(maize,maize_events,'pi',windows = 21)+theme(legend.position='none'),
plot_over_time(maize,maize_events,'xi',windows = 21)+theme(legend.position='none'),
plot_over_time(maize,maize_events,'tajD',windows = 21)+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/origVal_distance_summary_maize.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(maize,maize_events,maize_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/region_summary_maize.pdf',width = 12,heigh=9)

plot_ratios(maize,maize_events,windows = 21)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05

ggsave('figures/smooth_distance_summary_maize.pdf',width = 16,height=7)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.01265
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.00865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 7.4823e-05
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.01265
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.00865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 7.4823e-05

Generic models

Demographies

demog <- read.csv('demographies/generic_models.csv')
demog <- demog %>%
  gather(model,size,Model.1,Model.2,Model.3,Model.4,Model.5,Model.6,Model.7,Model.8,Model.9,Model.10,Model.11,Model.12) %>%
  mutate(bneck_length=case_when(model%in%c('Model.2','Model.3','Model.6','Model.7')~'1_Short',
                                model%in%c('Model.4','Model.5','Model.8','Model.9','Model.10','Model.11','Model.12')~'2_Long',
                                model=='Model.1'~'Constant'),
         bneck_strength=case_when(model%in%c('Model.2','Model.4','Model.6','Model.8')~'1_Weak',
                                  model%in%c('Model.3','Model.5','Model.7','Model.9')~'2_Strong',
                                  model%in%c('Model.10','Model.11','Model.12')~'3_constant',
                                  model=='Model.1'~'Constant'),
         bneck_age = case_when(model%in%c('Model.6','Model.7','Model.8','Model.9')~'2_Late',
                               model%in%c('Model.2','Model.3','Model.4','Model.5','Model.10','Model.11','Model.12')~'1_Early',
                               model=='Model.1'~'Constant'),
         size=size)
demog <- filter(demog,model!='Model.1')

length_labels <- c('1_Short' = 'Short','2_Long'='Long')
strength_labels <- c('1_Weak'= 'Weak','2_Strong'='Strong','3_constant'='Constant')
age_labels <- c('1_Early'='Old','2_Late'='Recent')


ggplot(demog,aes(generation,size,color=model))+ 
  geom_hline(yintercept = 20000,color='grey80')+
  geom_line(size=1.5) +
#  geom_hline(yintercept = 400)+
  scale_y_log10() +
  theme(axis.text.x = element_text(angle = 45,hjust = 1,vjust = 1,face = 'bold')) +
  labs(x=expression(paste('Generation (x ',N[anc],')')),
       y='Population size',
       color='Model')+
  facet_grid(bneck_age~bneck_length+bneck_strength,space = 'free',
             labeller = labeller(bneck_length=as_labeller(length_labels),
                                 bneck_age = as_labeller(age_labels),
                                 bneck_strength = as_labeller(strength_labels)))+
  theme(strip.background = element_blank(), strip.placement = "outside")
## Warning: Removed 72000 rows containing missing values (geom_path).

ggsave('figures/generic_models/generic_demographies.pdf',width = 12)
## Saving 12 x 5 in image
## Warning: Removed 72000 rows containing missing values (geom_path).

model1

model1_events <- c(0,0.9)
model1_event_names <- c(0,0.9)

Neutral data

model1_neutral <-read_sim("results/generic/model1_neutral_pi.csv",
                             "results/generic/model1_neutral_xi.csv",
                             "results/generic/model1_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model1_neutral,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model1_bgs <-read_sim("results/generic/model1_bgs_pi.csv",
                         "results/generic/model1_bgs_xi.csv",
                         "results/generic/model1_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model1_bgs,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model1 <- merge_sel_neutral(model1_bgs,model1_neutral)
head(model1)
##      gen stat window centimorgen      bgs  neutral  mean_val
## 1 -0.002   pi      1     0.00000 7.746545 13.31417 0.5818270
## 2 -0.002   pi     10     0.00738 8.398654 13.22437 0.6350891
## 3 -0.002   pi     11     0.00820 8.453766 13.27309 0.6369101
## 4 -0.002   pi     12     0.00902 8.593112 13.23504 0.6492697
## 5 -0.002   pi     13     0.00984 8.593202 13.17840 0.6520674
## 6 -0.002   pi     14     0.01066 8.640976 13.17171 0.6560254
var_pi <- function(theta,n){
  
  my_var <- (((n+1)*theta)/(3*(n-1)))+((2*(n^2+n+3)*theta^2)/(9*n*(n-1)))
  
}
model1_pi_var <- var_pi(13.28,400)
sqrt(model1_pi_var)/sqrt(1000)
## [1] 0.2093724
theta <- 13.28
n=400


neutral_pi <- read.csv("results/generic/model1_neutral_pi.csv")
head(neutral_pi)
##      gen    X0    X1    X2    X3    X4    X5     X6     X7     X8     X9
## 1 -0.008 8.130 6.455 7.739 4.571 6.413 9.515 11.278 23.822 27.641 22.982
## 2 -0.005 7.494 6.346 7.435 4.393 5.506 8.245 10.167 23.088 27.386 22.314
## 3 -0.002 7.480 6.145 7.494 4.478 5.670 8.304  9.948 22.789 26.858 22.266
## 4  0.000 7.707 5.784 6.982 4.600 6.162 8.587 10.491 23.800 27.238 22.731
## 5  0.002 7.785 6.102 7.222 4.618 6.089 8.676 10.636 24.006 27.406 23.177
## 6  0.005 7.442 6.119 7.363 4.299 5.653 8.100 10.081 23.655 27.590 23.674
##      X10    X11    X12    X13    X14    X15    X16    X17    X18    X19
## 1 18.017 15.643 17.159 19.201 32.965 24.574 18.818 12.661 10.248 12.846
## 2 18.403 15.911 18.317 19.211 33.172 24.907 20.046 13.443 10.486 13.373
## 3 18.228 16.044 17.707 19.227 32.763 24.874 20.149 13.378 10.799 13.164
## 4 18.238 15.721 18.088 19.238 32.783 24.502 19.329 13.527 10.915 13.006
## 5 18.014 15.779 17.278 19.180 32.910 24.607 18.657 13.024 10.842 12.909
## 6 18.195 15.877 17.559 18.757 32.426 24.228 19.491 13.388 10.862 13.215
##   replicate
## 1      2511
## 2      2511
## 3      2511
## 4      2511
## 5      2511
## 6      2511
tail(neutral_pi)
##           gen     X0     X1     X2     X3     X4    X5     X6     X7
## 2019995 0.988 17.257 11.415 11.729 10.725 10.443 8.615 13.411 17.169
## 2019996 0.990 17.136 11.437 11.878 10.838 10.095 7.859 12.341 15.283
## 2019997 0.992 17.798 11.555 11.609 10.771 10.732 9.078 14.224 17.401
## 2019998 0.995 18.433 12.233 12.231 11.150 10.783 8.834 14.301 17.977
## 2019999 0.998 18.240 11.904 11.882 10.980 11.125 8.966 14.447 18.557
## 2020000 1.000 18.560 12.118 12.088 11.218 11.242 9.199 14.752 19.128
##             X8     X9    X10    X11    X12    X13    X14    X15    X16
## 2019995 18.808 16.767 15.301 11.958 16.248 12.965 11.644 10.300 12.018
## 2019996 17.854 15.532 14.030 11.308 14.759 11.708 10.579  9.313 10.955
## 2019997 19.322 17.499 16.009 12.349 15.984 13.000 11.906 10.267 11.816
## 2019998 19.395 17.315 15.766 12.227 16.132 12.092 10.855  9.599 11.235
## 2019999 19.688 17.851 16.180 12.410 16.128 12.221 11.489  9.981 11.905
## 2020000 19.856 17.729 16.685 12.892 17.316 13.248 12.467 10.995 12.573
##            X17   X18   X19 replicate
## 2019995 11.199 8.521 6.381      4308
## 2019996 11.025 8.197 6.068      4308
## 2019997 12.624 9.186 6.088      4308
## 2019998 11.467 8.690 6.455      4308
## 2019999 11.533 8.802 6.414      4308
## 2020000 11.824 9.311 6.627      4308
neutral_pi %>%
  gather(value,window,-replicate,-gen) %>%
  group_by(replicate,value) %>%
  summarise(window=mean(window))%>%
  filter(window>13.28+sqrt(model1_pi_var)*2) 
## # A tibble: 3,060 x 3
## # Groups: replicate [1,460]
##    replicate value window
##        <int> <chr>  <dbl>
##  1         2 X3      31.6
##  2         2 X6      26.8
##  3         4 X10     27.5
##  4         4 X9      29.4
##  5        10 X10     28.0
##  6        10 X4      33.0
##  7        10 X5      29.8
##  8        10 X6      27.6
##  9        10 X7      34.4
## 10        13 X15     33.7
## # ... with 3,050 more rows
3060/100000
## [1] 0.0306
length(unique(neutral_pi$replicate))
## [1] 5000
ggplot(neutral_pi,aes(gen,X1,group=replicate))+geom_line()+
    geom_ribbon(aes(ymax = 13.28+sqrt(model1_pi_var)*2, ymin = 13.28-(sqrt(model1_pi_var)*2)))

neutral_pi %>%
  group_by(replicate)%>%
  summarise_at(c("X1", "X2"), mean, na.rm = TRUE) %>%
  filter(X1 >13.28+sqrt(model1_pi_var)*2, X2 >13.28+sqrt(model1_pi_var)*2) %>%
  count()
## # A tibble: 1 x 1
##       n
##   <int>
## 1    72
72/5000
## [1] 0.0144
mean(neutral_pi$window)
## Warning in mean.default(neutral_pi$window): argument is not numeric or
## logical: returning NA
## [1] NA

Plot absolute values

legend <- get_legend(plot_over_time(model1,model1_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model1,model1_events,'pi')+theme(legend.position='none')+geom_hline(yintercept = 13.28),
plot_over_time(model1,model1_events,'xi')+theme(legend.position='none')+geom_hline(yintercept = 13.28),
plot_over_time(model1,model1_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model1.pdf',width = 16,height = 12)

Check if neutral sims are stabile

ggplot(filter(model1,stat=='pi',window%in%c(1,10,20)),aes(gen,neutral,group=window))+
    geom_ribbon(aes(ymax = 13.28+sqrt(model1_pi_var), ymin = 13.28-sqrt(model1_pi_var)))+
  geom_ribbon(aes(ymax = neutral + sd_neutral/sqrt(1000), ymin = neutral - sd_neutral/sqrt(1000),fill=as.factor(as.numeric(window))),alpha=.2) +geom_line()+
  geom_hline(yintercept = 13.28)+

  labs(fill='window ')

ggsave('test/origVal_distance_summary_model1.pdf',width = 10,height = 10)

Plot ratios

simplot_grid(model1,model1_events,model1_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model1.pdf',width = 12,heigh=9)

plot_ratios(model1,model1_events)

ggsave('figures/generic_models/smooth_distance_summary_model1.pdf',width = 16,height=7)

model2

model2_events <- c(0,0.5)
model2_event_names <- c('Bottleneck','Recovered')

Neutral data

model2_neutral <-read_sim("results/generic/model2_neutral_pi.csv",
                             "results/generic/model2_neutral_xi.csv",
                             "results/generic/model2_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model2_neutral,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model2_bgs <-read_sim("results/generic/model2_bgs_pi.csv",
                         "results/generic/model2_bgs_xi.csv",
                         "results/generic/model2_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model2_bgs,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model2 <- merge_sel_neutral(model2_bgs,model2_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model2,model2_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model2,model2_events,'pi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'xi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model2.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model2,model2_events,model2_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model2.pdf',width = 12,heigh=9)

plot_ratios(model2,model2_events)

ggsave('figures/generic_models/smooth_distance_summary_model2.pdf',width = 16,height=7)

model3

model3_events <- c(0,0.6667)
model3_event_names <- c('Bottleneck','Recovered')

Neutral data

model3_neutral <-read_sim("results/generic/model3_neutral_pi.csv",
                             "results/generic/model3_neutral_xi.csv",
                             "results/generic/model3_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model3_neutral,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model3_bgs <-read_sim("results/generic/model3_bgs_pi.csv",
                         "results/generic/model3_bgs_xi.csv",
                         "results/generic/model3_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model3_bgs,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model3 <- merge_sel_neutral(model3_bgs,model3_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model3,model3_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model3,model3_events,'pi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'xi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model3.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model3,model3_events,model3_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model3.pdf',width = 12,heigh=9)

plot_ratios(model3,model3_events)

ggsave('figures/generic_models/smooth_distance_summary_model3.pdf',width = 16,height=7)

model4

model4_events <- c(0,0.05,0.5250)
model4_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model4_neutral <-read_sim("results/generic/model4_neutral_pi.csv",
                             "results/generic/model4_neutral_xi.csv",
                             "results/generic/model4_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model4_neutral,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model4_bgs <-read_sim("results/generic/model4_bgs_pi.csv",
                         "results/generic/model4_bgs_xi.csv",
                         "results/generic/model4_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model4_bgs,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model4 <- merge_sel_neutral(model4_bgs,model4_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model4,model4_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model4,model4_events,'pi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'xi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model4.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model4,model4_events,model4_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model4.pdf',width = 12,heigh=9)

plot_ratios(model4,model4_events)

ggsave('figures/generic_models/smooth_distance_summary_model4.pdf',width = 16,height=7)

model5

model5_events <- c(0,0.05,0.6833)
model5_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model5_neutral <-read_sim("results/generic/model5_neutral_pi.csv",
                             "results/generic/model5_neutral_xi.csv",
                             "results/generic/model5_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model5_neutral,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model5_bgs <-read_sim("results/generic/model5_bgs_pi.csv",
                         "results/generic/model5_bgs_xi.csv",
                         "results/generic/model5_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model5_bgs,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model5 <- merge_sel_neutral(model5_bgs,model5_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model5,model5_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model5,model5_events,'pi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'xi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model5.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model5,model5_events,model5_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model5.pdf',width = 12,heigh=9)

plot_ratios(model5,model5_events)

ggsave('figures/generic_models/smooth_distance_summary_model5.pdf',width = 16,height=7)

model6

model6_events <- c(0.,0.05)
model6_event_names <- c('Bottleneck','Recovered')

Neutral data

model6_neutral <-read_sim("results/generic/model6_neutral_pi.csv",
                             "results/generic/model6_neutral_xi.csv",
                             "results/generic/model6_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 220000
simplot_grid(model6_neutral,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model6_bgs <-read_sim("results/generic/model6_bgs_pi.csv",
                         "results/generic/model6_bgs_xi.csv",
                         "results/generic/model6_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 220000
simplot_grid(model6_bgs,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model6 <- merge_sel_neutral(model6_bgs,model6_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model6,model6_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model6,model6_events,'pi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'xi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model6.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model6,model6_events,model6_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model6.pdf',width = 12,heigh=9)

plot_ratios(model6,model6_events)

ggsave('figures/generic_models/smooth_distance_summary_model6.pdf',width = 16,height=7)

model7

model7_events <- c(0.0,0.0667)
model7_event_names <- c('Bottleneck','Recovered')

Neutral data

model7_neutral <-read_sim("results/generic/model7_neutral_pi.csv",
                             "results/generic/model7_neutral_xi.csv",
                             "results/generic/model7_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 220000
simplot_grid(model7_neutral,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model7_bgs <-read_sim("results/generic/model7_bgs_pi.csv",
                         "results/generic/model7_bgs_xi.csv",
                         "results/generic/model7_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 220000
simplot_grid(model7_bgs,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model7 <- merge_sel_neutral(model7_bgs,model7_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model7,model7_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model7,model7_events,'pi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'xi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model7.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model7,model7_events,model7_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model7.pdf',width = 12,heigh=9)

plot_ratios(model7,model7_events)

ggsave('figures/generic_models/smooth_distance_summary_model7.pdf',width = 16,height=7)

model8

model8_events <- c(0.0,0.05,0.0750)
model8_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model8_neutral <-read_sim("results/generic/model8_neutral_pi.csv",
                             "results/generic/model8_neutral_xi.csv",
                             "results/generic/model8_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 220000
simplot_grid(model8_neutral,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model8_bgs <-read_sim("results/generic/model8_bgs_pi.csv",
                         "results/generic/model8_bgs_xi.csv",
                         "results/generic/model8_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 220000
simplot_grid(model8_bgs,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model8 <- merge_sel_neutral(model8_bgs,model8_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model8,model8_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model8,model8_events,'pi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'xi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model8.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model8,model8_events,model8_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model8.pdf',width = 12,heigh=9)

plot_ratios(model8,model8_events)

ggsave('figures/generic_models/smooth_distance_summary_model8.pdf',width = 16,height=7)

model9

model9_events <- c(0.0,0.05,0.0833)
model9_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model9_neutral <-read_sim("results/generic/model9_neutral_pi.csv",
                             "results/generic/model9_neutral_xi.csv",
                             "results/generic/model9_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 220000
simplot_grid(model9_neutral,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model9_bgs <-read_sim("results/generic/model9_bgs_pi.csv",
                         "results/generic/model9_bgs_xi.csv",
                         "results/generic/model9_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 220000
simplot_grid(model9_bgs,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model9 <- merge_sel_neutral(model9_bgs,model9_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model9,model9_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model9,model9_events,'pi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'xi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model9.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model9,model9_events,model9_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model9.pdf',width = 12,heigh=9)

plot_ratios(model9,model9_events)

ggsave('figures/generic_models/smooth_distance_summary_model9.pdf',width = 16,height=7)

model10

model10_events <- c(0.0)
model10_event_names <- c('Bottleneck')

Neutral data

model10_neutral <-read_sim("results/generic/model10_neutral_pi.csv",
                             "results/generic/model10_neutral_xi.csv",
                             "results/generic/model10_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model10_neutral,model10_events,model10_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model10_bgs <-read_sim("results/generic/model10_bgs_pi.csv",
                         "results/generic/model10_bgs_xi.csv",
                         "results/generic/model10_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model10_bgs,model10_events,model10_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model10 <- merge_sel_neutral(model10_bgs,model10_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model10,model10_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model10,model10_events,'pi')+theme(legend.position='none'),
plot_over_time(model10,model10_events,'xi')+theme(legend.position='none'),
plot_over_time(model10,model10_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model10.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model10,model10_events,model10_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model10.pdf',width = 12,heigh=9)

plot_ratios(model10,model10_events)

ggsave('figures/generic_models/smooth_distance_summary_model10.pdf',width = 16,height=7)

model11

model11_events <- c(0.0)
model11_event_names <- c('Bottleneck')

Neutral data

model11_neutral <-read_sim("results/generic/model11_neutral_pi.csv",
                             "results/generic/model11_neutral_xi.csv",
                             "results/generic/model11_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model11_neutral,model11_events,model11_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model11_bgs <-read_sim("results/generic/model11_bgs_pi.csv",
                         "results/generic/model11_bgs_xi.csv",
                         "results/generic/model11_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model11_bgs,model11_events,model11_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model11 <- merge_sel_neutral(model11_bgs,model11_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model11,model11_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model11,model11_events,'pi')+theme(legend.position='none'),
plot_over_time(model11,model11_events,'xi')+theme(legend.position='none'),
plot_over_time(model11,model11_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model11.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model11,model11_events,model11_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model11.pdf',width = 12,heigh=9)

plot_ratios(model11,model11_events)

ggsave('figures/generic_models/smooth_distance_summary_model11.pdf',width = 16,height=7)

model12

model12_events <- c(0.0)
model12_event_names <- c('Growth')

Neutral data

model12_neutral <-read_sim("results/generic/model12_neutral_pi.csv",
                             "results/generic/model12_neutral_xi.csv",
                             "results/generic/model12_neutral_tajD.csv",
                             neutral=TRUE)
## [1] 2020000
simplot_grid(model12_neutral,model12_events,model12_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model12_bgs <-read_sim("results/generic/model12_bgs_pi.csv",
                         "results/generic/model12_bgs_xi.csv",
                         "results/generic/model12_bgs_tajD.csv",
                         neutral=FALSE)
## [1] 2020000
simplot_grid(model12_bgs,model12_events,model12_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model12 <- merge_sel_neutral(model12_bgs,model12_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model12,model12_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model12,model12_events,'pi')+theme(legend.position='none'),
plot_over_time(model12,model12_events,'xi')+theme(legend.position='none'),
plot_over_time(model12,model12_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model12.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model12,model12_events,model12_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model12.pdf',width = 12,heigh=9)

plot_ratios(model12,model12_events)

ggsave('figures/generic_models/smooth_distance_summary_model12.pdf',width = 16,height=7)